Effect of population size in a Prey-Predator model 

Fabien Campillo* Claude Lobry* 
November 29, 2011 



Abstract 



We consider a stochastic version of the basic predator-prey differen- 
tial equation model. The model, which contains a parameter uj which 
represents the number of individuals for one unit of prey - If x denotes 
the quantity of prey in the differential equation model x = 1 means 
that there are uj individuals in the discontinuous one - is derived from 
the classical birth and death process. It is shown by the mean of sim- 
ulations and explained by a mathematical analysis based on results 
in singular perturbation theory (the so called theory of Canards) that 
qualitative properties of the model like persistence or extinction are 
dramatically sensitive to uj. For instance, in our example, if uj = 10 7 
we have extinction and if uj = 10 8 we have persistence. This means 
that we must be very cautious when we use continuous variables in 
place of jump processes in dynamic population modeling even when 
we use stochastic differential equations in place of deterministic ones. 
Keywords: Prey-predator model; Ordinary Differential Equations; 
Diffusion Equations; Gillespie algorithm; Birth and Death processes 

1 Introduction 

Consider the standard prey-predator model : 



where x stands for the concentration of preys and y for the concentration 
of predators. It is well known that this kind of modeling with differential 
equations is valid only if one unity of x (or y) represents a large number 
of prey (or predator) individuals. On the other hand, when the number of 
individual is too small, everybody agree that one must switch to some kind 
of individually based modeling of stochastic nature. 
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What means large is generally not specified but it is widely admitted that 
around 10 3 the law of large numbers begins to do its job and that figures 
like 10 6 are completely safe if one wants to use continuous variables and 
differential equations. 

The objective of this paper is to show that the threshold of 10 3 is not always 
acceptable and that, in some circumstances, even 10 6 cannot be considered 
as secure when we deduce biological consequences, like persistence, from the 
behavior of a model with continuous variables. For that purpose we propose 
a stochastic model, where the dynamic of the prey is governed by a birth 
and death process while, for mathematical simplicity, we keep the predator 
variable as a continuous one. The development will make clear that this 
simplification does not affect the conclusions of the paper. The proposed 
model is such that the dynamic of the process is locally approximated (when 
the number of preys is large) by a differential system which is precisely of 
predator-prey type like 0. We agree that, in many respects, our model is 
biologically questionable but our objective is not to contribute to biologi- 
cal understanding of prey-predator relationship but just to point out some 
mathematical phenomenon which is likely to be present in many models and 
which might be responsible for erroneous interpretations. 

The first section is devoted to the presentation of the stochastic model, 
the second to the presentation of some surprising simulations, the third to 
the analysis of the differential system that governs the dynamics of the mean 
of the stochastic process and the forth to the explanations of the surprising 
aspects of the simulations. The last two sections are devoted to method- 
ological and bibliographical comments. 

From the mathematical point of view the material and results presented 
here are classical. The paper is intended principally for non mathemati- 
cally oriented readers who are not necessarily aware of these questions. We 
tried to avoid all mathematical technicalities and for this purpose we made 
an important use of results from computer simulations. All the references 
to existing literature related to these questions are rejected to the last two 
sections. 

2 The model. 

The variable uj x{t) is an integer which is the number of preys at time t. 
This variable performs the following birth and death (actually here "death" 
means "capture" by a predator) process. 

• At any time, the epoch r of the next event (birth or death) is a random 
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variable Z which follows an exponential law of parameter 



X = ^(f(x) +f x(x)y) (2) 

• At the epoch r we have one birth with probability f^^^ y or one 
death with the complementary probability : 

f{<r-)) 



P{ujx{t + ) =ujx(t~) + 1) 
P(cvx(t + ) =ujx{t~) - 1) 



f(x(r ))+/i(x(t ))y(r~ 
fx(x(r-))y(T-) 



(3) 



f{x{r-))+^x{r-))y{r-) 

The variable y is a continuous variable which evolves according to : 

y(t + dt) = y(t) — dtm y(t) + ^{number of captures during [£, t + dt]} (4) 

Thus the predator dynamics is an exponential decay associated to a growth 
proportional to the number of prey disappearing during the elapsed time. 
The parameter e accounts for different time scale for the prey and the preda- 
tor dynamics. 

Assume that dt = 10 -4 , uj = 10 9 , e = 10 _1 and f(x) + l*>{x)y is of the 
order of unity. Then, during elapsed time dt the number of events (death or 
birth) is of the order of A dt = j(f(x) + fi(x)y)dt w ^ w 10 6 . This is a bit 
lengthy to simulate (at least with a desk computer) but, due to that great 
number of events, the process defined by ([2]), Q, Q is accurately approxi- 
mated on the interval [t, t + dt] by the diffusion process (see appendix A for 
a derivation) : 

{ x(t + dt) = x(t) + f[f(x(t))-n(x(t))y(t)}-a x W t 

\ (5) 

[ y(t + dt) = y(t) + dt[(n(x(t)) - m)y{t)\ + a y W t 

where Wiat, W2M, W^dt, is a sequence of independent Gaussian variables 

with mean and standard deviation 1 with : 




f(x(t)y(x(t))y(t) 
f(x(t)) + v(x(t))y(t) 

f(x(t)n(x(t))y(t) 
f(x(t)) + v(x(t))y(t) 

This diffusion process is not a good approximation of the jump process when 
x is small. For an accurate description one must switch to the jump process 
for small values of x but, since this is not our point here, we restrict us to 
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the consideration of the stochastic (with continuous variables) diffusion like 
process : 

if x{t) < i then x(t + dt) = else 
x(t + dt) = x(t) + f[f(x(t))-fi(x(t))y(t)]+a x W t ( 6 ) 

y(t + dt) = y(t) + dt[(fi(x(t)) - m)y{t)} + a y W t 

The first line in ^ states that when the number of prey is smaller than 1 
it has to be . It is necessary to specify this because now the variable x is 
continuous in the diffusion process but we want to keep the meaning of x as 
a number of individuals . Thus, ^ must be an absorbing barrier for For 
x > ^ one sees that the recurrence equation for the mean of x(t) and y(t) 
is approximated by : 

f E[x(t + dt)} = E[x(t)] + f[f(x(t))-^x(t))y(t)} 

\ E[y(t + dt)} = E[y(t)} + dt[(n(x(t) - m)y{t)} 1 1 

which is the Euler scheme for the differential system : 

= \[f(x) - v( x )y] 

(8) 

Thus, to conclude this paragraph, we have constructed a diffusion-like model 
defined by equations Q. This model depends on a parameter uo. This model 
has the following properties : 




Since the model is derived from the birth and death process lux is 
interpreted as the number of individuals for x units of preys. 



• The standard deviation is proportional to y ^ : the biggest is uj the 
more "deterministic" is the process. 

• The diffusion process is degenerate (i.e. the dimension of the random 
noise is not 2 but 1). This is due to the fact that only x is considered 
as a discrete variable, not y. 

• When uj x is large (greater than 10 3 ) the dynamic of the model is 
accurately approximated (at least for small durations) by the classical 
deterministic differential prey-predator model ([§]) which is the same as 
system with c = £, em — 5 after a change of time units. 

We shall first simulate this system and then explain the observed simula- 
tions. 
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Figure 1: Above : uo = 10 12 ; below 
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3 Simulations 

In this section we fix /, /i, e and m as : 

• f( x ) = \ x ( 2 - x) 

• ji{x) 



0.4 + x 



• e = 0.02 



• rn = 0.6645 

On Fig |l] one sees one run of the process Q. The duration is 200 and uo 
is fixed at 10 12 (above) and 10 10 (below). One sees regular oscillations for 
the population of prey (in red) and predator (in black). We do not see any 
difference between the two records. These regular oscillations are those pre- 
dicted by the deterministic prey-predator model. Since the value of x during 
the oscillations is around 1 which corresponds to such a great number of in- 
dividuals we are definitely not surprised that the continuous deterministic 
system is a good approximation. 

But, on Figj2]we observe a dramatic change with uj — 10 8 which is still a big 
figure. We observe a mixed mode oscillation with a random successions of 
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Figure 2: Above : uj = 10 8 ; below uj = 10 6 

large and small oscillations which could not be produced by a deterministic 
two dimensional system. With uj — 10 6 we observe an extinction of the two 
populations which is confirmed on Figj3] where we observe that none of 20 
runs for uj — 10 6 is persistent at time T — 200. 

Let us denote by T the time of extinction for the predator (defined as the 
time where y(t) reaches the value ^). Let us say that there is extinction 
when the time of extinction is smaller than 1000. On Tab. Q]we have the 
empirical probabilities of extinction with respect to uj (computed on 1000 
runs) and mean and standard deviation of T computed on trajectories end- 
ing with extinction for t < 1000 . We can see that the transition is very 
sharp from extinction with probability one (uj — 4.0 10 6 ) to non extinc- 
tion with probability one (uj = 2.0 7 ). It seems surprising that with about 
uj = 2.0 10 7 prey-individuals the system is definitely (say up to 1000 units of 
time) safe and definitely unsafe for 4.0 10 6 which is still a big figure. This is 
a problem since in most case, in population dynamics models, we have poor 
information on the actual size of a population. We come back later on this 
issue. Notice also that the standard deviation of T is very large for small 
values of uj which makes predictions very imprecise. 
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200 f 

Figure 3: Twenty runs with uj — 10 6 

4 The dynamics of the continuous deterministic 
model 
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Figure 4: 

In this section we describe the dynamics of the deterministic model Q which 
approximate the evolution of the mean of the diffusion model Q. All the 
material in this section is classical and known as the theory of "canards" 
(see the section "literature comments" for more details). 

The first step in the understanding of a planar system like our is to draw 
the two nullclines (sometimes called "zero growth isoclines"), that is the sets 
defined by : 

• The nullcline of the prey : {(x,y) : ^ [/(#) — l^(x)y] = 0} 

• The nullcline of the predator :{(x,y) : (ii(x) — m)y = 0} 

In our simulations the parameter e is small (0.02) and, by the way, except 
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u 


E[T\ 


a(T) 


P(T < 1000) 


10 5 


30.46 


6.75 


1 


10 6 


39.02 


11.30 


1 


2.0 10 6 


47.74 


19.62 


1 


4.0 10 6 


79.05 


51.79 


1 


6.0 10 6 


143.54 


121.42 


0.999 


8.0 10 6 


259.76 


222.64 


0.983 


9.0 10 6 


311.70 


247.58 


0.964 


10 7 


554.21 


319.94 


0.867 


1.1 10 7 


555.17 


351.75 


0.741 


1.2 10 7 


681.31 


324.12 


0.649 


1.3 10 7 


745.83 


321.95 


0.481 


1.4 10 7 


815.46 


296.26 


0.384 


1.5 10 7 


867.54 


273.60 


0.255 


1.6 10 7 


906.10 


238.55 


0.182 


1.7 10 7 


928.68 


221.50 


0.120 


1.8 10 7 


964.05 


143.82 


0.072 


1.9 10 7 


975.48 


110.09 


0.059 


2.0 10 7 


> 1000 








Table 1: Empirical probabilities of extinction according to uo 
when the quantity 

[f(x) - fi(x)y] 

is small, of the order of £, the right member in the first equation in Q is 
large compared to the second one. This means that the vector velocity of 
([§]) is almost horizontal. From this it follows that, a first approximation the 
solutions of our system is shown by the hand-drawn schemes on Figj4] and 
Figj5] : Outside of the parabola and the y axe which is the nullcline of the 
prey the trajectories are taken as horizontal. 

• On Figj4] one sees that the nullcline of the predator (in blue) is on the 
left of the maximum of the nullcline of the prey (the black parabola 
curve). Along the nullcline of the prey the motion is down-up on the 
right of the vertical blue nullcline and up-down on the left. From this 
we see that there is a tendency for the trajectories to join the y-axe 
on its attractive part (above the black curve), to follow it in the up- 
down direction and there is some indeterminacy to where it will leave 
it after having crossed the the parabola. From this scheme we suspect 
the existence of a periodic limit-cycle cycle which, actually, can be 
proven to be the case. 

• On Figj5] the situation is somewhat easier to understand. The blue 
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Figure 5: Schematic representation of solutions of ^ 

vertical null-cline of the predator being on the right the motion along 
the parabola converges to a limit point which apparently is a stable 
attracting point for all initial conditions. 

• Notice that an attracting equilibrium and an attracting limit-cycle are 
qualitatively different picture and that the transition between the two 
cases occurs when m crosses the value 0.666.... (when the blue line 
crosses the parabola at its maximum). 





Figure 6: Phase portrait of system Q for different values of m 

Let us now comment on Figj6} The pictures are not hand-drawn schemes 
but actual simulations with e = 0.02 ; we observe the great similarity with 
the schemes. 
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• rn = 0.6 : One large limit cycle (the direction of the motion is counter 
clockwise). Trajectories above the limit cycle are of two kinds : some 
hit the limit cycle and then follow it and the others hit the vertical axe, 
then they follow it up-down and reappear below, run left to right hit 
the parabola and then join the limit cycle. Actually "true" trajectories 
never meet but, due to the limit of our drawing, they seem to meet. 
All trajectories follow for a while the y axe and then x(t) is potentially 
small. 

• m = 0.75 : We have an attracting equilibrium. Some trajectories go 
directly to the equilibrium, some other follow the y axe. 

• m — 0.6645 : We have a small periodic limit cycle circling around the 
unstable equilibrium which is very close to the periodic orbit. Along 
the periodic orbit we stay very far from the y axe (and, by the way, 
x(t) is never small) but one sees that near the unstable equilibrium, a 
very small perturbation leads to a trajectory which follows the y axe 
and x{t) can become small. 

• m — 0.66442561 : In this case we have a limit cycle which is of inter- 
mediate size between "large" (follows the y axe for a while and small 
(remains far fro the y axe) ; it just hits the y-axe. The point is that it 
needs very sharp values for m (8 digits in our case) to obtain this inter- 
mediate cycle called a "canard cycle" . See in the section "comments" 
some informations about the mathematical theory of "canards" . 

All along this description we said that x(t) is potentially small when the 
trajectory follows the y axe. But how small ? A simple way to enlarge what 
is going on along this axe is to plot, not the point (x(t), y(t)), but (£(£), y(t)) 
with : 

C(t) = e\n(x(t)) 

This is done on Figj7|and Figj8j We represent the (x, y) and and the (£, y) 
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Figure 7: (x, y) and (£, y) variables on the same axes : m = 0.6 
trajectories in the same system of axes ; (x, y) trajectories are in red, (£, y) 
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Figure 8: (x,y) and (£, y) variables on the same axes : m = 0.6445 

are in green and both limit cycles in the two systems of representation in 
blue. The two vertical red lines correspond to x = 10 -9 and 10 -6 . There 
are 19 trajectories starting from (2,0.5 ± A; 0.05) k = 0, 1, ...,9. 

Let us compare the two simulations. 

• Fig|7| We look at the "large" limit cycle in the (£, y) variables and we 
see that the minimum of £ corresponds to x — 10 -9 ; for the trajectory 
labeled 1 the minimum is about 10 -17 . These incredibly small values 
are easily explained in Appendix B. 

• Figj8j The "small" limit cycle is almost not visible in the (£, y) vari- 
ables. In both the (£, y) and the (x, y) variables the trajectories labeled 
1 to 4 look very similar. But, in Figj8j trajectory 5 remains above 10 -9 
which is not the case in FigjTjand trajectory 6 does not exist in FigjT] 

The main difference between the case m — 0.6 and the case m — 0.6645 
is that, in the first case, every trajectory is such that the minimum of x 
is smaller than 10 -9 unlike in the second case where there are two set of 
trajectories : Those that start above trajectory 6 for which the minimum 
will be smaller than 10 -6 before reaching the limit cycle and the others for 
which x{t) remains greater than 10 -6 . Notice that this trajectory 6 is in 
some places very close to the limit cycle. 

The observed differences between m = 0.6 and m = 0.6645 are not specific 
of these values. In particular the same behavior with two type of trajec- 
tories separated by a sharp transition is true for all values of m between 

m — 0.66442561 and m — 0.6666 This behavior is summarized by the 

description of the "safety funnel" shown in Figj9] by the green arrow and 
that we explain now. Assume that for some reason we do not accept to 
pursue a trajectory such that the min of x(t) is smaller than a = 10~ k (it 
may be because we think that the size of the population is to low in order 
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to survive or because we want to switch to a different - stochastic - model). 
The form 10~ k is by no mean essential for a, it is just to emphasize that a 
is small. It exists a unique yo such that the solution issued from (oo, y$) (in 
practice 2 is a good infinite), which we call the "a-safety trajectory", is such 
that x(t) first decreases and attains a first local minimum equal to 10 _fc . 
This is the red trajectory on the scheme of Figj9j This trajectory, when 
x{t) « ^ _1 (m) will be very close to the limit cycle (the blue trajectory) . 
We call p(e, k) the distance between the two curves ; this can be evaluated 
from the value of e and k. The "safety funnel" is defined by the parts of 
red and blue curves on the right of the vertical x = /x _1 (m). If a trajectory 
which enters the funnel is perturbed, as long as it remains in the funnel, 
the (future) minimum of x will remain greater the 10~ k . If not, there is a 
danger to reach values smaller than 10~ k . 
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Figure 9: The "safety funnel' 



5 The diffusion process in the variables (x,y) and 



UJ 


10 9 


10 8 


10 7 


10 6 


p(e,k) 


1.2 10" 3 


9.0 10" 5 


5.5 10" 5 


5.3 10" 5 


a x \fdt 


4.2 10" 5 


1.4 icr 5 


4.2 10" 5 


1.4 10" 4 


(j y \fdt 


4.9 10" 9 


1.4 10" 7 


4.9 icr 7 


1.4 10" 6 



Table 2: Width of the funnel and corresponding a x and a, 
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On the four simulations shown on FigfT0|to 13 we have performed 10 runs of 
20 time units duration of the process ^ starting from (—2, 0.5). The results 
are presented in both (x,y) and variables (black trajectories). In the 
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same variables we have simulated from system ([8]) the "safety trajectory" 
corresponding to 1000 individuals, that is to say : The "safety trajectory" 
was obtained by dichotomy and the observed width p(e, k) of the funnel 
is given on the table [2] with corresponding rough evaluations of a x and a y 
around the funnel. 

• FigjlO} All the runs are widely below the "safety trajectory". 

• Fig jlll All the runs are below the "safety trajectory" but we observe 
that some runs are close to it. 

• Fig|l2j All the runs are above the "safety trajectory" and some are 
closse to the vertical line corresponding to 10 individual. 

• Fig. |12| All the runs are widely above the "safety trajectory" and 
reach ultimately the vertical line corresponding to 1 individual. 

We observe that when uo decreases the strength of the randomness increases 
and at the same time the width of the funnel decreases. These opposite 
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Figure 12: cj = 10 7 




Figure 13: u = 10 6 

trends are responsible for the sharp transition from extinction to persistence 
as uj grows from 4.0 10 6 to 2.0 10 7 . 



6 Methodological comments 

6.1 About the question of size of populations. 

We are used to the fact that continuous differential models works rather well 
in fluid dynamics and chemical kinetics despite the ultimate discrete nature 
of fluids. We know that this efficiency is related to the very large number of 
atoms in the process. Von Foerster, Lotka, Volterra and others popularized 
the formalism of chemical kinetics in the domain of population dynamics 
; they were certainly aware of the limits of such an approach but, in the 
absence of computers and with a far less developed probability theory, it 
was a way to progress. 



Now, thanks to computers and probability theory, we have good models 
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for small populations. Unfortunately these models are still expensive in 
term of computer time and deterministic or diffusion models (stochastic 
differential equations) with continuous variable are still unavoidable. In a 
diffusion model the size of the population considered is directly related to 
the strength (standard deviation) of the random term. 

The example of prey-predator interaction presented here shows that the qual- 
itative behavior of such models may depend strongly on the size of the pop- 
ulation even when it is very large. 

6.2 About the generality of the example 

The deterministic prey-predator model Q approximate the dynamics for 
E[x] and E[y] of the birth and death model defined by ([2]), ([3]), Q and its 
diffusion approximation Q. This model Q is the very classical determin- 
istic prey-predator model which is proposed in every text book as a first 
improvement of the Lotka-Volterra model. The separation of time scales for 
prey and predator dynamics introduced by the presence of the parameter e 
in the first equation has the following classical explanation. Using a change 
of time unit (Tsl) rewrites : 



If we use the same mass unit for x and y then e is a yield factor. A yield 
factor like 0.02 is acceptable in ecology (one needs 50 kg. of dry grass to get 
1 kg. of cow). For bigger s like 0.1 the sharp transition that we presented 
is still present but less spectacular. 

As previously said we admit that our birth and death model is question- 
able with respect to its biological signification. There are certainly many 
different models for individual behavior with the same deterministic equa- 
tion approximating the mean of the process. Since our point relies on the 
diffusion approximation for such models our conclusions are valid as long as 
such approximation is correct. In the case of birth and death processes it 
works provided that the number of individuals is greater than 10 3 -10 4 which 
is our case. For more elaborated models at the individual scale (for instance 
physiologically structured preys) this point remains to be considered. 






em — 8 
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6.3 About the existence of "canard " solutions in the model. 



Let us say two words about "canard solutions " . In a system with two time 
scales like : 



consider the curve T defined by the equation y) = ; this curve split in 
two regions : 

• the attracting one made of points such that, in the neighborhood, the 
vector field converges to T, 

• the repelling one made of points where, in the neighborhood, the vector 
field diverges from T, 

separated by equilibria. A "canard " solution is a solution of the differen- 
tial system which follows, for some duration, the attracting part of T at a 
distance of the order of e and, after that, follows also the repelling part at 
a distance of order s. Some " canards " are robust which means that they 
persist under small changes in the model, others are not. 

The presence of a "safety funnel" like the one described in section 3 is 
related to the presence of two "canard solutions" in Q. 

• The solution t — >> (x(t) = 0, y(t) = y(0)e~ mt ) which corresponds to the 
absence of prey, 

• a solution following the cubic from the right to the left, which has no 
analytic expression but which existence can be proved by continuity 
arguments. 

The first "canard" is robust but the second is not. This is the reason why, 
the sharp transition between 4.0 10 6 and 2.0 10 7 individuals occurs for a 
rather short interval of values of the parameter m. As a consequence, to 
some extend, our example is exceptional, not "generic". This will be the 
case in most two dimensional systems, but this do not invalid our point 
since robust "canard" (different from trivial "canard" corresponding to the 
absence of some population) are generically present for dimension 3 and 
more. An easy way to understand it is to imagine that our parameter m is 
of the form : 



which mimics, for instance, some seasonal dependence of the mortality rate. 
This non autonomous system can be considered as a three dimensional sys- 
tem and we see that the "canard" value for m is crossed periodically. We 




m(t) = a + b cos(r t) 
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Figure 14: u = 10 12 10 9 (above left, right) uj = 10 8 10 7 (below left, right) 

have done a simulation in the case : 

m(t) = 0.6645 - 0.047(1 - cos(O.H)) 

and the results are shown on Fig[l4j For uj = 10 12 we observe no difference 
between the deterministic model {x in blue, y in black) and the diffusion 
approximation (x in red, y in black) ; for uj — 10 9 we observe a very slight 
deviation between red and blue curves ; for uj — 10 8 we observe a very big 
difference with now a mixed mode oscillation in the diffusion process ; for 
uj — 10 7 the mixed mode oscillation leads to extinction. 

6.4 About the inadequacy of deterministic models with con- 
tinuous variables. 

In population dynamics every body agrees that deterministic models are just 
crude approximations of reality. Only individually based models, stochastic 
by essence, can represent correctly the evolution of real ecosystems. The 
example presented here is just one more argument against the danger of 
using deterministic differential equations without care. 

But it is by no mean an argument against the study of continuous deter- 
ministic differential models of populations dynamics ! 

Actually there are many good reasons for continuing to explore systems 
of ordinary differential equations : 

• Some models are mathematically appealing. For instance the proof 
of the exclusion principle for the most general model of competition 
in the chemostat [12], despite its poor ecological contents, remains an 
interesting mathematical challenge for mathematicians. 
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• More interesting is the use of easily tractable mathematical models to 
formalize some ecological issue and clarify the discussion. An inter- 
esting example of this use of differential equations is given by the dis- 
cussion on "ratio dependent" models initiated by the paper of Ardity 
and Ginsburg pQ. 

• In our example, the understanding of the diffusion model, relies on 
very particular and recently (see bibliographical comments) discov- 
ered properties of deterministic differential systems : the "canard so- 
lutions" . 

By the way, far from being an article of propaganda against the use of 
deterministic differential systems, our paper supports the importance of a 
thorough understanding of the properties of ordinary differential systems in 
population dynamics. In particular it shows that the classical deterministic 
definition of persistence : 

limsupx(t) = a > 
must be enriched by some consideration about the "size" of a. 

6.5 About computer simulations in dynamic population mod- 
eling. 

There is no doubt that our mathematical understanding of the phenomena 
outlined in the present paper will considerably increase in the future. But 
this will require high mathematical sophistication and time. Unfortunately, 
in the mean time, biologist will use models and computer simulations which 
are not completely safe. It urges to provide them with computer routines 
which are safe of numerical artifacts associated to the true nature of a popu- 
lation : a more or less large number of individuals. Considering our present 
mathematical knowledge this certainly can be done in a comparatively short 
time but it needs quite a lot of people working on the design of safe com- 
puter software. This was done in the past for the needs of industry (for 
instance digital wind tunnels), medicin (medical imaging) this could be the 
case for microbial ecology but it depends of decisions at the level of scientific 
policies. 

7 Bibliographical comments. 
7.1 The atto-fox problem. 

The question of the inadequacy of deterministic continuous modeling is 
firmly addressed by D. Mollison [10J in a paper which criticize the biological 
interpretations of a previous paper by Murray et al. [TT] . Let us quote from 
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As to the second wave, close inspection shows that the expla- 
nation lies, not much in the determinism of the model, as in its 
modeling of the population as continuous rather than discrete 
and its associated inability to let the population variables reach 
the value zero. Thus the density of infected at the place of ori- 
gin of the epidemic never becomes zero, it only declines to a 
minimum of around one atto-fox (10 -18 of a fox, Hugues 1960) 
per square kilometer. The model then allows this atto-fox to 
start the second wave as soon as the susceptible population has 
regrown sufficiently. 

About ten years before Mollison, independently, within the framework of 
chemical kinetics, D. Gillespie published a famous paper |6j : Exact Stochas- 
tic Simulations of Coupled Differential Reactions from which our model in 
the present paper is inspired. 

It is a bit surprising that, at least to our knowledge, not much has been 
done in this direction. The present paper is a development of a first draft [9J 
with T. Sari where we noticed the importance of the presence of "canards 
solutions" regarding the question of persistence in ecological models. The 
paper [7J is also related to this atto-fox question in the case of the chemostat 
with a slow varying flow rate. The paper [3] which is much more mathe- 
matically oriented, considers the stochastic modeling of the chemostat ; it 
focusses on the the approximation of jump processes by diffusion processes 
and was a source of inspiration for the present paper. 

7.2 Singular perturbations and "canard solutions". 

As already said, "canards" are specific solutions in singular perturbations of 
differential equations. They where discovered in 1981 by a group students of 
G. Reeb : E. Benoit, J-L. Callot, F. and M. Diener [2j. They studied them 
within the framework of Non Standard Analysis which is most suitable for 
modeling since it is a simple formal language where the use of infinitesi- 
mals (in the sense where physicists use this term) is mathematically rigor- 
ous. But they are now also studied by numerous mathematicians within 
the framework of mached asymptotic expansion or the geometric singular 
perturbation theory. The article [13J by Martin Wechselberger is a short and 
nice introduction to "canards" and the paper [5] is a thorough survey about 
our present understanding of "canards" with a focus on numerical questions. 
The paper [8] is about Nonstandard Analysis applied to real word questions. 

The question of considering the presence of noise in singularly perturbed 
systems has been considered for long time. We refer to the recent paper [3J 
devoted to the question of the consequence of noisy environment on "canard 
solutions" and its bibliography. In particular the results contained in this 
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paper allow to give asymptotic evaluations of the wide of the "safety funnel" 
and many other quantity of interest but their mathematical sophistication 
is out of the scope of the present paper. 



8 Conclusion 

Scientists are now much familiar with the phenomenon of "sensitivity to 
initial conditions" which, in some deterministic dynamical systems, is the 
cause of an impredictable long range behavior. The same phenomenon in 
some deterministic differential equations modeling the dynamic behavior of 
populations is the cause that a very small difference in an initial condition (or 
along a trajectory) will make the future value of some variable very small 
or not. This is the reason why, in the modeling of population dynamic, 
it is a good thing to add some small noise to the deterministic process 
because it does not cost too much computer time and may detect this kind 
of phenomenon. But we have shown that the result may depend strongly 
of the strength of the noise. By the way, when we do not have an accurate 
estimation of the strength of the noise, it should be more secure to vary that 
strength and make sure that the behavior is not strongly dependent on it. 



9 Appendix 

A Approximation by a diffusion process 

Consider the process defined by ([2]), ([3]), ([I]). Since Z follows an exponential 
law of parameter A its expectation is^ and the number iV& of events during 
the duration dt is approximately : 

N b w dt/{\) =dt\ = dt-(f(x)+fjL(x)y) 

A £ 

We consider iV& as deterministic. If dt is small the variables x(t) and y(t) 
are approximatively constant. Denote by AQ the random variable which is 
equal to one if at the i-th event a predation occurs ; one has : 

P( Xi = +l) Kx(t))y(t) 



p(x = o) = fix(t)) 

[ 1 ' f(x(t))+n(x(t))y(t) 

The number of predations during [£, t + dt] is, approximately, J2^ b AQ and 
the number of birth is by the way N b — ^2i b AQ and the increment of the 
number of individuals is N b — 2 J2^ b X{. 



One has : 
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Figure 15: Comparison between birth and death process and diffusion pro- 
cess 



.E[X<]- 



f(x(t)) + ii(x(t))y(t) 



E[J2X 1 } = dt^(f(x(t)) + rtx(t))y(t)) " {xit)) y[t) 



^ £ v,v v „ ^ v n f( x (t)) + n(x(t))y(t) 

dt^ n(x(t))y(t) 
m 2, Y , f(x(t))fi(x(t))y(t) 



{f{x{t))+^x{t))y{t)Y 
o\ J> ) = dt^f(x(t)) + ,(x(t))y(t)) fW))Mt))v{t) 



j e w ^-" ' — — (/(x(t)) + /i(x(t))y(t))2 



s (f(x(t)) + fi(x(t))y(t)) 

From the central limit theorem we can approximate the sum by a Gaussian 
and we write : 



V e ^wm; y £ ( / ( x ( t )) + /Li ( ;r ( t ))y( t )) 

where Wi is a Gaussian of mean and 1 as standard deviation . 

Since the variable x is the number of individuals divided by uo the incre- 
ment of x is given by : 



x(t + dt) -x{t) = -(N b -2 



and replacing by the value of Nf, one get : 
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x{ t+dt)-x{t) * dt hf(x(t))-n(x(t))y(t)}-Jdt- f/f^f^ t $£m Wi 

e y cue (/(x(t)) +M»(*))f(*)) 

Let us compute now the increment of y. According to Q we have : 

y(t + dt) — y(t) — dtm y(t) + ^{number of prey death during [t, t + dt] } 
which, according to the previous notations is : 

y(t + dt) - y(t) = -dt m y(t) + ~Y j X i 
and introducing Wt one gets : 

y( t + *) - yit) « W ) ~ nWQ + ft J^f^S )™' 

On Fig|l5] one sees a comparison between the birth and death process (red 
trajectories) and its approximation by a diffusion. From the left to the right 
we have uj — 10 6 , uj — 10 5 , uj — 10 4 . The representation is both in (x,y) 
and variables. We have 10 runs from the initial condition (0.2, 0.6). 

The two red vertical lines correspond to a population between 1 and 1000. 



B Exponentially small values 

Let us write explicitly system Q as : 
dx 1 



dt 
dy 
dt 



[0.5x(2 -x)- — 

e v v ; 0.4 + x 



y] 



x 



v 0.4 + x 

In the variables (£, y) the system writes : 



m)y 



1 



r ft = ^ 2 -^-oi+m 



y] 



dy 
dt 



(; 



my 



0.4 + i/e 

which is approximated, when £ << e, by : 

ft - 

dy 

Tt = ~ my 



(10) 



(ii) 
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Take as initial condition (xq, yo) = (—0.1, 0.9) (which corresponds to trajec- 
tory n° 1 in Fig. [7| and integrate. It comes that the minimum for ^{t) is 
attained for the value t* of t for which y(t*) = 0.4 and this value turns out 
to be approximately —1. But : 

x * = e & « e - 40 « 10- 17 

The minimum depends much of the value of yo : The largest is yo the smallest 
is the minimum. This explain why in FigjH] the minimum corresponding to 
trajectory 6 is much bigger. 

C Numerical simulations 

We did not use any solver. A specific software was written in order to be 
sure that there were not artifacts caused by erroneous uses of some sophis- 
ticated numerical scheme. Trajectories of the differential equations ([§]) are 
obtained using the Euler scheme defined by Q. We prefer this scheme to 
any more sophisticated scheme used to simulate differential systems since it 
is the exact recurrence scheme which approximate for E[x(t)] and E[y(t)] of 
the diffusion process Q. 

We fixed dt = 10 -4 since we observed that for this value solutions of Q 
are indistinguishable from those with dt = 10 -5 . 

The birth and death process defined by ([3]), Q takes too long time 
to be simulated when A is very large (in the case of our computer uj > 10 6 ) 
and this is the reason why we used a diffusion approximation which is a 
perfect approximation for large values. Since we where mainly interested 
in the funnel phenomenon associated to "canard " it was not necessary to 
switch to the true birth and death process for small values of A. But if one is 
interested by figures like the mean of the extinction time it should be better 
to switch to some suitable jump process. 
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